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Abstract 



The dynamics of a 2D site percolation model on a square lattice is studied us- 
ing the hierarchical approach introduced by Gabrielov et al., Phys. Rev. E, 60, 
5293-5300, 1999. The key elements of the approach are the tree representation of 
clusters and their coalescence, and the Horton-Strahler scheme for cluster ranking. 
Accordingly, the evolution of percolation model is considered as a hierarchical in- 
verse cascade of cluster aggregation. A three-exponent time-dependent scaling for 
the cluster rank distribution is derived using the Tokunaga branching constraint 
and classical results on percolation in terms of cluster masses. Deviations from 
the pure scaling are described. An empirical constraint on the dynamics of a rank 
population is reported based on numerical simulations. 

1 Introduction 

Percolation model is probably the simplest and best studied system that experiences 
(geometrical) phase transition of the second kind 34j. It is widely used as a toy model 
for spatially distributed stochastic processes, such as diffusion in disordered media, for- 
est fires, gelation, semiconduction, etc. jSHISll- Importantly for our study, percolation 
model presents a transparent mechanism of the process of hierarchical aggregation (coag- 
ulation). This process has been actively employed for describing the essential properties 
of material fracture and earthquake nucleation [U IHl 1101 lEl EH 123 EHl EHj , starting from 
the pioneering works of AUegre et al. [Tj and Newman and Knopoff [221 12ni ^] • In this 
paper we describe the evolution of percolation model in terms of an inverse cascade of 
hierarchical cluster aggregation. 

An early idea of hierarchical aggregation was introduced by Newman and Knopoff in 
the "crack-fusion" model for repetitive cycles of large earthquakes [221 123 El Ell • Their 
model focused on processes of small cracks fusions into successively larger ones, accommo- 
dating the influence of mainshocks and aftershocks, juvenile crack genesis from tectonic 
stresses, crack healing, and anelastic-creep induced time delays, plus other effects. Tur- 
cotte et al. have reinstated this line of research considering a log-binned description 
of hierarchical aggregation and performing numerical tests to study its scaling proper- 
ties. Gabrielov et al. ^T] first have employed the Horton-Strahler hierarchical ranking 
[ini I2n| to construct an exactly solvable model of a general inverse cascade process. The 
Horton-Strahler ranks (see Sect. 12. 3|) that came from hydrology and have been not well 
known in physical applications happened to be more natural than cluster masses (sizes, 
areas) in describing the aggregation process. Moreover, the ranks are shown essential for 
formulating the analytical models [TT] . Recent efforts deal with studying the aggregation 
dynamics and its various scalings via exactly solvable hierarchical models and extensive 
simulations |18j . 

Below we focus on the evolution of the first spanning cluster in the the classical site- 
percolation model, and decribe it as a consecutive hierarchical fusion of smaller clusters 
into larger ones. Noteworthy, we are interested not in a final solution of a percolation 
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state, but in an evolutionary path leading from the juvenile single-particle clusters to 
a self-similar population of clusters of arbitrary large size (limited by the finiteness of 
the lattice), the percolation cluster included. Thus we depart from the steady-state 
assumption of [HI UHl EH] as well as from the asymptotic focus on the percolation onset 
typical for the classical percolation studies |31j. 

Specifically, we follow [TT] and represent each cluster by a time-oriented tree that 
reflects the history of cluster formation. The model dynamics is then described in terms 
of the corresponding trees using the well-developed theory of hierarchical scaling com- 
plexities 121 EE]- An important role is played by the the Horton-Strahler scheme that 
provides a natural ranking for the tree-based structures. Another important element 
is the Tokunaga classification that defines a special subclass of trees with self-similar 
branching. A large number of hierarchies observed in nature are shown to be Tokunaga 
trees j26] ; this is also the case for the clusters in percolation model |^ . We use the 
Tokunaga constraint together with classical results on percolation dynamics (in terms 
of cluster masses) to derive time-dependent scaling laws for rank distribution of clus- 
ters. Importantly, we report a three-exponent scaling for the dynamics of a population 
of clusters of a given rank, in deviation from the two-exponent scaling well-known for the 
population of a given mass |5H ITT] . We also analyze deviations from the pure scaling 
and confirm our results by numerical simulations. 

The inverse cascades and aggregation (coagulation) processes are important for evo- 
lution of many natural hazardous processes: earthquakes, landslides, and forest fires are 
argued to follow the hierarchical aggregation dynamics [101 ^] . A general review of the 
theory and models of kinetics of irreversible aggregation is given by Leyvraz ^Hl- An 
alternative approach to analytical modeling, based on ideas from IT], but using equa- 
tions that are consistent with the mass action law of chemical kinetics, can be found in 
da Costa et al. [7]. 

The paper is organized as follows. The percolation model is described in Sect.|21 this 
section also introduces tree representation of clusters and the Horton-Strahler ranking. 
In Sect. 01 we derive the average mass of clusters of a given rank using the Tokunaga 
constraint on cluster branching. This result will be actively used in consecutive sec- 
tions. Section E) is devoted to the time-dependent rank distribution of clusters. First 
(Sect.mH), we establish the exponential rank distribution at percolation using the result 
of Sect. El We then proceed with time-dependent rank distribution; Sect. 14.21 introduces 
the three-exponent scaling for ranks and compares it to the well-known Stauffer's two- 
exponent scaling for cluster masses. Scaling for ranks averaged over the entire evolution 
of the percolation cluster is derived in Sect. 14.31 this result is motivated by the heuristic 
studies that typically use averaged observations on a system. Time-dependent finite-size 
corrections to the established scalings are described in Sect. 14.41 Our study of rank 
distributions is concluded in Sect. 14.51 bv describig the time-dependent behavior of the 
total mass of clusters of a given rank. Sect. El analyzes fractal properties of clusters and 
reports sharp increase of cluster fractal dimension in the vicinity of percolation. Sect. El 
uses simulations to establish a notable constraint on the dynamics of a rank populations. 
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The results are discussed in Sect, d 



2 Model 



2.1 Dynamics 

We consider the classical 2D site-percolation model [SI]. The model dynamics starts with 
an empty Lx L square lattice. At each step a particle is dropped into a randomly chosen 
unoccupied site; thus each site can be either occupied by one and only one particle or 
empty. Two sites are considered neighbors if they share one side; each site on a square 
lattice has four neighbors. Cluster is defined as a group of occupied neighbor sites p^ . 
Time refers to the steps at which particles drop onto the lattice. Since we do not have 
annihilation of particles, time is formally equivalent to the number of particles on the 
lattice. It is convenient to normalize time by the lattice size so it varies from p = 
at the start to p = 1 when all sites are occupied. During the system evolution, occupied 
sites start to aggregate and clusters begin to form. Once a sufficient number of particles 
is accumulated, a percolation cluster is formed connecting the opposite sides of the lattice 
vertically and/or horizontally. 

The density p increases monotonically from zero to its critical value pc at percolation. 
For an infinite lattice pc ~ 0.59274606 |25j, while for a finite lattice it is smaller jMj : 



Many phenomena encountered in the percolation model mimic what we see when the 
phase transitions of the second kind occur. Note however that these phenomena are 
of purely geometrical and statistical rather than physical nature. Indeed, the physical 
percolation theory is largerly predicated in this geometrical model and there are many 
empirical links between them; this is why the percolation model is said to be an example 
of the geometrical phase transition of the second kind, and why its nomenclature emerges 
from that of the physical critical phenomena. 

The theoretical description of the percolation dynamics is conventionally given in 
terms of the cluster masses [S^! and most of the universal scalings - a benchmark of phase 
transitions of second kind - deal with parameters expressed via the mass distribution 
of clusters. However, if one is interested in analytical description of the aggregation 
process, the mass description happens to be inferior to the hierarchical rank approach 
fU^]. Properly defined ranks not only allow one to construct exactly solvable models 
of aggregation, but also they are more feasible for observations in practice. In addition, 
they refiect the individual history of cluster formation. Below we follow the hierarchical 
approach of Gabrielov et al. to study the percolation dynamics. 



Pc{L) = pc-cL 
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2.2 Tree representation of clusters 

Each cluster in our model is represented by a tree that reflects the time- dependent for- 
mation of a cluster (its history), and is a subject for quantitative analysis. Specifically, 
each one-particle cluster is represented by a trivial tree consisting of a single node. When 
two clusters are merged together their trees are also merged by adding a new node (par- 
ent) for which they become children (and siblings to each other.) In our model, the 
coalescence of two or more clusters can only be materialized by adding to the lattice a 
new particle which will be a neighbor to one or more existing clusters. Figure ^ illus- 
trates the four possible types of coalescence. We call fc-coalescence in a situation when a 
newly dropped particle (marked N in the figure) is a neighbor to k existing clusters (gray 
numbered sites). Numerical simulations on a square lattice with L = 2,000 suggest the 
following relative frequencies Qk of fc-coalescences: Qi ~ 0.628, Q2 ~ 0.318, Q3 ~ 0.052, 
Q4 ~ 0.002. Figures ^,c illustrate how a tree is formed for different coalescence types. 
There are two basic situations: When a new particle is a neighbor to only one existing 
cluster, it is considered as an individual one-particle cluster that is connected to the 
existing one. The connecting node of the tree in this case does not correspond to a 
particle on the lattice (panel b). When a new particle is droped in a neighbor to two, 
three, or four existing clusters, it is not condidered as an individual cluster. Instead, it 
corresponds to the connecting node in the tree (panel c). Thus, the connecting node in a 
tree may or may not correspond to a lattice particle depending on the coalescence type. 
The branching parameter (number of children for a given parent) of a tree for any cluster 
varies between 2 and 4. Note that both 1- and 2-coalescences result in merging only two 
clusters; accordingly, most of the observed coalescences (about 95%) involve only two 
clusters while coalescence of three or four clusters is extremely rare. 

The consecutive process of tree formation for a simple four-particle cluster is illus- 
trated in Fig. 121 Importantly, the individual evolution of a cluster is crucial in construct- 
ing the corresponding hierarchical tree. To construct the tree one needs to consider all 
consecutive coalescences that have formed the cluster, not only its final shape. Therefore, 
it is clear that the same tree may correspond to clusters of different shape: Figure Et' 
shows two 11-particle clusters that both correspond to the same tree shown in panel b. 
Therefore, working with trees, we unavoidably narrow the information about the cluster 
population. Notice however both trees capture an excessively larger amount of informa- 
tion than mere cluster masses. Summing up, the time evolution of a cluster is neccesary 
and sufficient to uniquely determine the corresponding tree, while the inverse is not true. 
The problems of describing the set of trees that might correspond to a given cluster, 
and the set of clusters that correspond to a given tree is beyond the scope of this paper. 
Next, we describe the ranking of clusters, presenting a conventional alternative to the 
logarithmic binning of cluster masses. 
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2.3 Horton-Strahler ranking 



The appropriate ordering of trees (clusters) is very important for meaningful description 
and analysis of the model dynamics. The problem of such an ordering is not trivial 
since the clusters may grow and coalesce in a variety of peculiar ways. An advantageous 
way to solve this problem is given by the Horton-Strahler topological classification of 
ramified patterns [13 ESI I2| illustrated in Fig. Eb: One assigns ranks to the nodes of a 
tree, starting from r = 1 at leaves (clusters consisting of one particle.) When two or 
more clusters with ranks r^, i = 1, . . . ,n merge together, a new cluster is formed with 
the rank 




ri + 1, if rj = ri V i = 1, . . . , n 
max(rj), otherwise. 



The rank of a cluster is that of the root of the corresponding tree. It is possible to consider 
an alternative definition of ranks: When at least two clusters with rank r coalesce, and 
other coalescing clusters have a lower rank, the rank of a new cluster becomes r + 1. 
Clearly, the two definition coincide when only two clusters coalesce. The results reported 
in this paper are independent of the particular definition, since coalescence of more than 
two clusters (especially of high ranks) is a rare event. 

Originally introduced in geomorphology by Horton jT3] and later refined by Strahler 
[SS], this classification is shown to be inherent in various geophysical, biological, and 
computational applications CH dl I2H1 EH iH ■ 



3 Mass-rank distribution 

In this section we derive the distribution of the average mass m,. of rank r clusters. It 
will be used consequently to connect various mass and rank scaling laws. First, we define 
the branching ratio Tij for a given cluster (tree) as the number Nij of subclusters (nodes) 
of rank i that joined subcluster (node) of rank j, averaged over subclusters (nodes) of 
rank J [201 EE]: 

- AT, • 

Next we note that the mass of a rank r cluster is the sum of two r — 1 cluster masses 
that formed the cluster (we ignore the possibility for three or more clusters to coalesce at 
the same step), plus a unit mass of a joining particle, plus the mass of all the lower-rank 
clusters that joined the considered cluster, hence: 

mi = 1 

m2 = (2 mi + P) + ri2(mi + P) 

ms = (2m2 + l)+T23(m2 + l) + Ti3(mi + P) 

fc-i 

TTik = (2mfc_i + l) + ^Tfc_,,(mfc_, + l)-(l-P)Tifc, A;>3. (2) 

i=l 
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Here the coefficient P addresses the possibihty for a one-particle cluster to join another 
cluster in two ways: via a one-particle connector (with probability P) or directly (with 
probability 1 — P); the clusters with r > 2 can only join other clusters using a one-particle 
connector. 

It was predicted by Gabrielov et al. fl] and later con&med by simulations JH] that 
clusters in percolation model obey the Tokunaga scaling j^nj asymptotically in k: 



SqS 



k-1 



(3) 



This rewrites Eq. for A; > 3 as 



k-l 



ruk = (2 rrik-i + 1) + E Ti{mk-i + 1) - (1 - P)Ty 



1=1 



Assuming the mass-rank relation in the form = d , c > 1 we obtain 



fe-i 



.fc-i 



2c''-' + l + Y^ sos'-' {c'-'-' + l) - (1 - P)sos'-' 



i=l 



k-2 



,k-l 



1 1 — (s/c)" " Sn S 



s/c 



Qk 2 s — 1 



l-P)so{s/c 



.k-2 



It is easily checked that this equation has a solution only if c > s; thus s/c < 1 and for 
large k then follows 



2 + 



So 



1 - s/c 



leading to the final equation 



- c(2 + s + So) + 2s = 



with solution: 



2 + s + So ± V(2 + s + So)2 - 8s 



(4) 



Remarkably, the model of Gabrielov et al. ^1] predicts in a Euclidean (assuming 
clusters of regular, non-fractal, shape) limit of an inverse cascade model 

So ^ 0.55495813, s = 1/so ^ 1.80193774, and c = 1/sg ^ 3.24697602. 

The Eq. (0)) in this case gives c(so, s) = 3.24697960 (this is the only solution such that 
c > s), which is remarkably close (6 digits) to the result of ^Jjj. Furthermore, the non- 
Eucledian (assuming fractal shape of clusters) steady-state simulations of Morein et al. 
[IH] suggest 

s ^ 3.0253, So ^ 0.6993, c ^ 4.325, 
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which exactly solves Eq. Q. We found it quite amazing that our complimentary set of 
assumptions used to derive (H)) lead to the same numerical results as analytical study [TT] 
and simulations of ^Hj- This suggests an underlying connection between our approaches 
to describe the hierarchical aggregation. 

The observed mass-rank distribution of clusters at percolation is shown in Fig. |31 it 
obeys the exponential relation 

= lOT('^-i) = c'^-\ (5) 

with 7 ^ 0.625, c = 10"^ ~ 4.2. Our simulation suggest that the mass distribution within 
a given rank is approximately lognormal (not shown) with the mean given by Eq. (0) 
and a rank-independent standard deviation. 

The relation (jS]) is a key element in our further analysis. As we will show, the distri- 
bution of cluster ranks at percolation (Sect.HIH) and its finite-size corrections (Sect. 
are obtained from the corresponding classical laws for masses by simple substituting the 
relation (jSj). At the same time, one of the most important results: the time dependent 
rank distribution can not be obtained this way and requires an additional treatment 
(Sect.Q. 

The exponential relation of Eq. (0) happens to be valid over the entire time interval 
< p < Pc- The corresponding dynamics of c(p) is shown in Fig. it grows with 
time from about 2.0 at the earliest stages to 4.2 at percolation. This growth reflects the 
fact that clusters become more weighty with time due to coupling with the clusters of 
lower ranks (which does not change the rank but increases the mass) . The growth is not 
monotonous; it is accompanied by pronounced log-periodic oscillations which are associ- 
ated with creation of new ranks. The log-periodic oscillations that accompany general 
power-law increase of observed parameters have been found in many systems including 
hierarchical models of defect development [2Z], biased diffusion on random lattices 
diffusion-limited aggregation (DLA) and others. Log-periodic oscillations can be 
naturally explained by the Discrete Scale Invariance (DSI) [21], which occurs in a system 
whose observables scale only for a discrete set of values. A famous example of DSI is 
given by the Cantor set that pocesses a discrete scale symmetry: In order to superim- 
pose its scaled image onto the original, one has to stretch it by the discrete factors 3", 
= 1,2, . . ., not a continuous set of values. The Cantor set and percolation belong to 
systems with built-in geometrical hierarchy, leading to DSI. In our particular system, 
ranks take only a countable set of values. Creation of new ranks necessarily disrupt the 
system in a discontinuous way resulting in the log-periodicity. 

Now we return to the numerical value of parameter c. In steady-state simulations of 
[TH] c = 4.325, which is reasonably close to what we observe at percolation. Recall that 
the models of jTH HH] use the "fractal correction" e to the cluster shape; this correction 
affects the rate r^j of clusters coalescence: 
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where Li is the total boundary size of the clusters of rank i. The correction e can be 
expressed as 

1 c-1 

which, together with results of Fig. shows that in the percolation model e decreases 
in time passing the Euclidean limit e = 1 j2] at {pc — p) ~ 0.14 and approaching the 
steady-state "fractal" e ~ 0.68 [18j at p = pc- The interval 2 < c < 4.2 observed during 
< p < Pc corresponds to 0.68 < e < oo. 



4 Rank distribution 

This section is devoted to establishing various time-dependent scaling laws for clusters of 
a given rank. We will see that it is tipically impossible to derive such laws by applying 
the mass-rank relation (jSj) to the coresponding well-known laws for cluster masses. This 
illustrates an original character and richness of the rank description and prompts for de- 
veloping new methods of analysis. We start with the simplest problem: rank distribution 
at percolation. 

4.1 Distribution at percolation 

We start recalling the well-known cluster mass distribution at percolation [33]: 

Umipc) ^ qom"", (6) 

where UmiPc) is the number of clusters of mass m per lattice site, and the Fisher exponent 
r = 187/91 ~ 2.05 is universal for 2D systems jHlIM]- Figure 10 illustrates the mass distri- 
bution at percolation for a system with L = 2000; to smooth out statistical fluctuations 
it shows the number of clusters with mass equal to or larger than m: Y,m'>m'^m{Pc)- 
Equation © suggests the slope r — 1 ^ 1.05, while the observed slope 0.96 is somewhat 
less than that. This is due to the impact of two concurrent phenomena: so-called "de- 
viation from scaling" at small m and finite-size effects at large m pTj I14j: they are 
discussed below in Sect. 14.41 

Now, we use Eq. ^ to derive the distribution of the number rir^Pc) of the clusters 
of rank r at percolation. Taking summation over all clusters of rank r and mass m we 
obtain: 

mup 



T - 1 

go 

r- 1 



rrir / \ m. 



m;^+\ (7) 
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Our simulations suggest (not shown) that the mass distribution within a given rank is 
lognormal with a rank-independent standard deviation. Thus, for arbitrary upper and 
lower quantiles m-up, mio of this distribution the values 

^lo(up) 

rrir 

are rank independent. Using this, we finally express nr{pc) via rrir'- 

ririPc) = Po m;^^^ QC m;^-°^. (8) 

The power law ^ is observed in a steady-state aggregation model of [TH] with index 
1.147. This index increase comparing to our 1.05 is due to the fact that in jTH] interme- 
diate clusters are removed from the lattice providing extra space for a larger number of 
smaller clusters. 

Combining the mass-rank relation (0) with (jS)) we obtain the following exponential 
rank distribution at percolation: 

UriPc) ~ pom;^+'=po{c-^+'y=pilO-'' (9) 

with 

Pi = Po c^~\ 6 = (r - 1) logio c ~ 0.62. 

This is indeed what we observe in Fig. [7| where the rank distribution Ur at percola- 
tion is shown by the dash-dotted line. The study ^H] suggests c^~'^ = 0.186 while our 
predictions and observations lead to c^~^ ~ ^_2~^-o^ = 0.22. The two values are in good 
agreement, the slight difference is explained, as in Eq. (jHJ, by removal of intermediate 
clusters in [TH]. Next we consider the rank distribution for p ^ pc. 

4.2 Dynamical rank distribution: three-exponent scaling 

Here we expand results of the previous section by establishing the time-dependent rank 
dustribution. First, we consider the dynamics of rank population. 

4.2.1 Temporal dynamics of rank population 

The dynamics of the total number (n^ ■ L^) of the clusters of a given rank r is illustrated 
in Fig. ISlfor r = 5,6,7. The population follows a characteristic bell-shaped trajectory, 
with percolation at its rightward limb. As in the case of mass description, one does 
not observe steady-state behavior in the cluster dynamics: The population of each rank 
steadily develops to its peak as a result of merging of the clusters of lower ranks; then it 
starts decreasing, giving birth to the clusters of higher ranks. As naturally follows from 
the model definition, the peak of the population of a higher rank comes after the peak 
of a lower rank. Figure IHl shows the population dynamics for the ranks 1 < r < 11 in 
semilogarithmic scale. Here one clearly sees the similarity in the dynamics of different 
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ranks. Note that this figure is remarkably similar to Fig. 7 from [3^] that shows the 
dynamics of clusters with logarithmically binned masses. We now proceed by establishing 
the appropriate time- dependent scaling laws. 



4.2.2 Time-dependent mass distribution 

Recall that the temporal dynamics of the cluster mass distribution is given by the two- 
exponent scaling law EH Hlj : 

Umip) ~ fo{z), z = {pc- p)nf + Zq, (10) 

with a = 1/2. The function /o has a bell-shaped form with maximum to the left of 
percolation; it can be roughly approximated by a Gaussian function ^3 E! : 

/o(z) oc exp {-az^^ . (11) 

Note that the shift zq is independent of m. 

Considered as a function of m, the two-exponent scaling explains the power law mass 
distribution ® at percolation (with go = /o(^o)) as well as the downward bend for 
p < Pc, clearly observed in Fig. (dashed line); while as a function of p it describes the 
bell-shaped dynamics of clusters with given mass m. 



4.2.3 Time-dependent rank distribution 

Combining the scaling laws ^ and (jlOj) one formally obtains the two-exponent scaling 
for rank dynamics. However, the two exponent scaling does not work for ranks; to show 
this we assume more generally 

ririp) ~ 9o{z)lO-'\ ^ = (Pc - p)h{r) + z'^, (12) 

which is consistent with the exponential rank distribution of Eq.Q at percolation (with 
Po = goi^'o))- Possible deviations from the pure exponential law at p < pc (clearly 
observed in Fig. [7j) and dynamics of a given rank (see Figs. I8l9p are described by specific 
form of the functions go{-) and h{-). Following we define 

:= ^ = (13) 

and choose h{-) in such a way that positions of the peaks of Urlz) coincide for different 
r; it is always possible by choosing the appropriate time change h{r). Figure ITUk shows 
the ratio Uri^z) / uii^z) for r = 2,3,6,8. One can see that the two-exponent scaling does 
not work in our case: the curves do not coincide. 

Nevertheless, the simple scaling picture is restored by introducing the additional, 
third, shift exponent: 

/i(r) = ailO"^", z;(r) = aslO-"^". (14) 
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Function go still can be approximated by a Gaussian function 

go{z) oc exp f - ^ ) . 




(15) 



Once the correct scaling form is established, the use of (0) is again legitimate, and 
the exponent cxi in Eq. ()14|) can be evaluated as: 



where c ~ 3 is the median of c values observed during p < pc- The observed exponent 
0"! ~ 0.23 (not shown) is fairly close to its predicted value. The shift exponent is estimated 
as (72 ~ 0.03; while scale coefficients are ai ~ 1.54, 02 ~ 1.43. The function gQ{z) that 
uses these estimates is shown in Fig.lTTlwhere different symbols depict clusters of different 
ranks. The collapse is obvious, confirming the validity of the three-exponent scaling p2|) . 



In the scaling for cluster masses, the time renormalization (p^ — p)m'' collapses the 
dynamics of mass m clusters onto the master curve /o(-2 — zq) with its only peak shifted 
by Zq leftward from percolation; the shift Zq is mass independent. Similarly, in the scaling 
for ranks the time renormalization (pc — p)10'^^'' collapses the dynamics of rank r clusters 
onto the master curve go{z — Zq), although the shift now is rank dependent and is given 
by 10'^^''. To illustrate this, we show the position of percolation on the righthand limb 
of the Gaussian go{z — 0.51) in Fig. llOb. The higher the rank, the closer the position of 
percolation to the peak of go- 

4.3 Averaged scaling 

In applications, it is often impossible to measure the size distribution of system elements 
at a given time instant. Moreover, sometimes the instantaneous size distribution does 
not exist at all: This is indeed the case for the systems described by marked point 
processes widely used to model seismicity, volcano activity, starquakes, etc. 0. In such 
situations one uses the averaged measurements. For instance, the famed Gutenberg- 
Richter law [121 EH 0] that gives exponential approximation to the size distribution of 
earthquakes (via their magnitudes) is valid only after appropriate averaging over a wide 
spatio-temporal domain. This explains the importance of the question: How do the 
distributions of Eq. (fTUj). (fT^ change after temporal averaging? 

We answer this question for averaging over < p < pc- For the mass distribution this 
leads to: 



(Ti = (TlogioC ^ 0.24, 



(HI, (Hni). 



n, 



•m 



OC 



OC 



oc 




(16) 
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Here the last step neglects the weak dependence of the integral on m (and uses the 
values T ^ 2.0, a = 1/2). The validity of p6p is confirmed by the observed averaged 
mass distribution shown by the solid line in Fig. IHl The averaged mass distribution is 
similar to that at percolation: it retains the power-law form while the slope is increased 
by 1/2 due to averaging. 

Similarly, we obtain for ranks: 

:= / nj.{p)dp = / go{z) 10 ^dp 
Jo Jo 

rpc 



oc / exp^-a' (pe-p)lO"^" -alO-"^"J^jlO-^"cip 

Jui ^ ^ 

oc 10^'^'^^+^)'" = I0'''('^^'^^^)^°sioc _ 10"^"'-. (17) 

The exponent may vary from 0.71 to 0.93 depending on 3.0 < c < 4.2 (the range of c 
values for the time when at least three ranks have been formed so the estimation of the 
distribution slope is meaningful). Simulations suggest (solid hne in Fig. Ej) = 0.87, 
which is in good agreement with our prediction. Again, the averaged rank distribu- 
tion retains the exponential form of the distribution at percolation; while its index has 
increased due to averaging. 



4.4 Correction to simple scaling 

Due to finiteness of the lattice, the results of previous sections require some corrections 
to match exactly the simulated rank distributions. The appropriate corrections are de- 
scribed below. 



4.4.1 Corrected scaling at percolation 

The pure power and exponential laws in Figs. El 13 are just first-order approximations to 
the observed cluster distributions at percolation. In both cases one sees the downward 
bending for small clusters and upward bending for large clusters. These are not due 
to statistical fluctuations. The downward bending for small clusters is explained by 
"deviations from scaling" |T3]: it can be shown analytically that the small clusters do 
not yet obey the general scaling law of Eqs. (jHl), (jHl) which holds only for large enough 
masses (ranks). The upward bend at large clusters is due to finite-size effects [TH IT7j: 
each large cluster that reaches outside the lattice boundary is "seen" as a number of 
smaller clusters, thus creating the upward deviation from the pure power (exponential) 
law. This phenomenon is especially important when the system is close to percolation 
and clusters of arbitrary large sizes have already been formed. The appropriate scale 
corrections for the mass distribution were studied by Hoshen et al. and Margolina 
et al. [nj. 
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To study the above phenomena it is convenient to consider the normahzed functions 

m'>m 

which, in the absence of scale corrections, would become constants: 

r — 1 

The function is shown in Fig. IT^ : it clearly deviates from the horizontal plateau at 
both sides. 

In case of the mass distribution, the corrections to scaling are given by ^7] : 

nm{pc) ^ m-^ (go + qi mT^ + qL m^^^L-^) , (18) 

where Q ~0.75, 1/D = 48/91 is the universal mean cluster radius exponent, and qo,qi,qL 
are independent of s and L. The first additional term describes the deviation from scaling 
for small clusters, while the second one is responsible for finite-size effects. 

For rank distribution, the "deviations from scaling" at lower clusters are only observed 
for r = 1; while the finite-size effects at large clusters are clearly present for many ranks. 
Accordingly, we propose the following correction to scaling for the rank distribution: 

Uripc) ^ 10'''' {po + PL 10''' L-') , r > 1. (19) 

with 

d = ^ log,, c ^0.33. 

The observed value of d can be estimated by plotting {rir 10^*" — po) as a function of r as 
shown in Fig. 112b. The observed ranks 4 < r < 9 follow the predicted scaling ()19|1 nicely. 

Importantly, the corrections to scaling ()19|) act at all cluster sizes, so they can not be 
neglected even for the intermediate clusters, not only for the largest ones. Indeed, their 
effect decreases with L, but this decrease is very slow. Notably, as shown by Morein et 
al. ^Hj (their Fig. 5) even for lattices as large as L = 30, 000 during the process when 
clusters as large as 2% of the lattice size are removed, the cluster size distribution clearly 
exhibits the upward deviations at large ranks (r = 11, 12, 13.) For smaller systems these 
deviations become dominant and may lead to an artificial decrease of the observed slope 
of cluster size distribution; this is demonstrated in Fig. l6|7l and is also seen in the analysis 
of Turcotte et al. [211 (their Fig. 9). 

4.4.2 Dynamics of scaling corrections 

Since the finite size effects play an important role in shaping the observed cluster size 
distribution, it is worth studying their dynamics. Specifically, we will be interested in 
transition of the cluster size distribution from the convex shape (in semi- or bilogarithmic 
scale) at p -C pc to formation of the upward bend at percolation. 
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For this we introduce a measure of convexity for the rank distribution, defined as an 
area between log nr{p) and a chord connecting its first and last points as shown in Fig. 
(the point r = 1 is not considered being affected by the deviations from scahng): 

/"^max 

M := [logion,(p) - (Ar + B)] dr, (20) 

with 

A = — , B = logio ^2 - 2A. 

The values of n are positive when nr{p) is convex in semilogarithmic scale, negative when 
it is concave, and vanish when it is linear. The measure /i(p) averaged over 1,000 runs 
on the lattice L = 2000 is shown in Fig. El the bell-shaped form of fi is decorated by 
the logperiodic oscillations for (pc — p) > 10~^ explained by creation of new ranks, which 
temporarily increases convexity. Zero level is crossed at about (pc — p) = 2 ■ 10~^, after 
that the rank distribution is concave. A detailed analysis (not shown) demonstrates that 
the distribution is never exactly linear; the transition from convex to concave shape is 
realized through the wave-shaped form when the distribution is still convex for the lower 
r, but is already concave for the higher ones. Qualitatively the same picture is observed 
for the mass distribution nm{p) (in bilogarithmic scale). 

The transformation of the cluster size distribution prior to percolation is not unlike a 
well-known pattern "upward bend" first described by Narkunskaya and Shnirman [T^ l^ 
in an early static model of defect development. Later it was found in steel samples and 
seismicity of California [28J, and confirmed by the dynamical modeling of failure in a 
hierarchical system (so-called colliding cascade models) [TUlli^ . 

4.5 Mass dynamics of a given rank 

Here we consider the dynamics of total and average mass of rank r clusters: 

r = }_^mnrm, = — ^ = — . (21) 

Here rirm denotes the number of clusters of rank r and mass m. Figure IT^ shows n^, M^, 
and nir for rank 5; the similar picture is observed for other ranks. It is tempting to use 
Gaussian approximation for Mj. and predict the Gaussian dynamics of (as a ratio of 
two Gaussians) and relate their parameters. Detailed analysis however demonstrates that 
under this approach the peak of m,, for ranks r > 9 should be observed after percolation; 
while in simulations this peak is always prior to percolation (not shown). Note that one 
still might approximate and by Gaussians with properly scaled parameters; such 
approximations will be good for rough curve fitting, but will fail to reproduce deeper 
properties of cluster dynamics. This demonstrates the general limitations of Gaussian 
approximations in the percolation problem. 
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5 Cluster fractal structure 



In this section we evaluate the fractal structure of clusters considering the mass-circumference 
relation 

m oc (22) 

where C is the number of empty neighbors of a cluster of mass m. For percolation 
cluster at infinite grid we have D^, = 1? which shows that the percolation cluster is a 
"linear" rather than a space-filling object [34j. Figure ITHk shows the cluster masses as a 
function of their circumference for different ranks. It is easily seen how the linear scaling 
Dr = 1 gradually develops as rank increases. Figure ITHb shows the index Dr estimated 
for 1 < r < 9. 

Figure IT^ shows the dynamics of D^ prior to percolation; noteworthy, its steady state 
behavior is altered by a gradual increase as p ^ pc- Similar increase is observed for 
clusters of other ranks. 

To explain the increase of Dr we recall that the rate of cluster coalescence is propor- 
tional to their circumference (see e.g. ^I]). Thus, for a given mass, clusters with a lower 
Dr have larger circumference, and a higher chance to coalesce. When a sufficient number 
of rank r clusters have been formed, the clusters with low Dr start to coalesce leaving 
the high-Dr clusters on the grid. 

Another reason for the increase of D^ is the finite-size effects. Specifically, this is an 
effect of having clusters that on an infinite grid have already gained higher ranks, but on 
our finite lattice are still small. 



6 Dynamical constraint 

Here we report an interesting regularity in rank dynamics that put a notable constraint 
on analytical modeling of percolation process. Specificaly, we consider the slope between 
two consecutive points of the rank distribution: 

er{p) := log 



n^+i(p)' 



Dynamics of ^4 is shown in Fig. [T7k together with that of n^. Noteworthy, the peaks of 
two curves (minimum of ^4 and maximum of rig) coincide. This happens to be true for all 
ranks: positions of corresponding peaks are shown as a function of rank in Fig. 117b. Such 
perfect matching is very unlikely to be accidental. Thus we conjecture that in order for 
nr{p) to properly describe the time-dependent behavior of rank population, the following 
system of differential equations must have a solution: 

_^ I rir ?^r+l ~ rirfir+i = (23) 
I nr+2 = 
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Applying this constraint to the three-exponent scahng of Eqs. (fT ^ .(fT ^ .(fT3 j) we find 



According to (j21I), the observed value ai = 0.23 gives (J2 = 0.04, which is 33% larger 
than the observed value (T2 = 0.03. The discrepancy is due to the approximate character 
of the Gaussian approximation (jl5|) for qq. 

7 Discussion 

The goal of this study was to describe the evolution of percolation model in terms of 
consecutive aggregation of smaller clusters into larger ones using the Horton-Strahler 
hierarchical scheme. First, this contributes to a novel understanding of the percolation 
phenomenon as a time-dependent hierarchical inverse cascade process. Second, this allows 
one to test the validity of the approach introduced by Gabrielov et al. and further 
developed by Morein et al. [IH] for a steady-state approximation to a general aggregation 
process. 

We considered dynamical and scaling properties of site-percolation on a 2D square 
lattice. Following [TT] we described clusters by hierarchical trees that reflect the history 
of cluster formation; the Horton-Strahler scheme was used to rank the trees and thus 
the corresponding clusters. We concentrated on the development of the first percolation 
cluster, thus working with a system that does not exhibit the steady-state dynamics, con- 
trary to the studies [TTl ITH] that have developed mean-field steady-state approximations 
to the system. 

Combining the results obtained in the classical percolation studies with the Tokunaga 
constraint on the cluster branching structure we derived various rank-dependent scaling 
laws connecting the number rir of clusters of rank r, their average mass rrir, and the rank 
r. We have compared the parameters of these laws with those predicted and observed by 
im ITH] in steady-state aggregation models. The values of parameters are shown to be in 
a perfect agreement, confirming the validity of the approach used in [TH I18j. In absence 
of the steady-state behavior, we derived the time-dependent versions of the scaling laws. 
We reported the three- index scaling (fT^. (fT^ for the number nr{p) of clusters of a given 
rank, which deviates from the classical two-exponent scaling for masses. 

We studied in detail the transition of the system from earlier stages to the vicinity 
of percolation and reported several characteristic phenomena observed as p — * pc- They 
include transformation of the cluster size distribution not unlike that observed in seis- 
micity, steel samples, and previous models of hierarchical fractures jlHl EHl CHI ESj ; and 
increase of the cluster fractal dimension. In our simple model these phenomena are partly 
explained (qualitatively as well as quantitatively) by finite-size effects; nevertheless we 
believe that they should not be neglected as irrelevant side-effects of numerical simula- 
tion. In fact, in practice we often work with systems that are described by intermediate 
depth hierarchies (in other words they have intermediate number of degrees of freedom). 




(24) 
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The percolation results related to small and intermediate lattices might be of high rele- 
vance in describing such systems. In addition, simulations on large lattices (L = 30, 000) 
performed by Morein et al. show that finite-size effects are still present even for large 
L. 

We have formulated an empirical constraint of Eq. for the time-dependent be- 
havior of rank population size nr{p)] the constraint follows very clearly from the observed 
values of nj.{p). It would be interesting to check this condition in real systems tradition- 
ally described by the percolation model. 

Our closing remark is on the index r of cluster mass distribution at percolation 
(Eq. (jni)). The studies of Gabrielov et al. [TTj and Morein et al. ^] predict r = 2; 
which slightly deviates from the well established theoretical value of the Fisher exponent 
r = 187/91 ~ 2.05. The index of the mass distribution is an essential characteristic of a 
system, thus even this slight difference of 2.5% might seem disappointing implying the 
intrinsically approximate character of the modeling of fTJ^j. In fact, this implication 
is not true. To validate the approach of jTT| ^] we notice that the Fisher exponent is 
tightly connected to the precise count of cluster particles on a site-level, hardly feasible 
in practice. At the same time, the studies |^ |H] have demonstrated that when we 
"characterize the size distribution of clusters in a way that circumvents the site-level 
description" considering any "macroscopic measure of the length scale of the cluster", 
the exponent of the corresponding scaling law becomes 2, universally for all 2D systems. 
An example of a "macroscopic measure" is the linear size in arbitrary direction, the radius 
of gyration, the diameter of the covering disk, etc. Clearly, the modeling of [HllIHl deals 
with such a macroscopic measure of cluster size, and hence predicts the correct scaling 
exponent. 
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Qi^0.63 Q2^0.32 Qg^O.OS Q4<0.01 




Figure 1: Multiple coalescence of clusters, a) Coalescence of clusters is materialized by 
adding to the lattice a new particle N (black) that is a neighbor to one, two, three, or 
four existing clusters (numbered gray sites). The relative frequencies Qk , k = 1,2,3,4 
of /c-coalescences based on similations with L = 2, 000 are shown in the figure. The 
corresponding tree is constructed as shown in panel b) (for k = 1) and c) (for k = 2). 
The cases = 3,4 are analogous to = 2. Note that about 95% of coalescences result 
in merging two clusters. See text for details. 
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Figure 2: Tree representation of clusters: scheme. The dynamics is from left to right. At 
first step particle A is dropped onto the lattice and a one-particle cluster is formed; it is 
represented by a one-node tree. At second step another one-particle cluster B is formed; 
it is represented by another one-node tree. At third step new particle C coalesces with 
cluster A to form two-particle cluster AC. This cluster is represented by a three-node 
tree; note that the connecting node of the tree does not correspond to any particle. At 
fourth step new particle D connects existing clusters AC and B to form four-particle 
cluster ABCD. This cluster correspond to a five-node tree. 



a b 




Figure 3: a) Non- uniqueness of tree representation. Two different 11-particle clusters that 
correspond to the same tree shown in panel b). Particles have been dropped according 
to their alphabet marks; so first was the particle A, then B, etc. b) Horton-Strahler 
ranking: illustration. The ranks are shown next to the tree nodes. 
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Figure 4: Mass-rank distribution observed on a 2,000x2,000 lattice at percolation. Dots 
- individual clusters, balls - average mass within a given rank. Line shows the relation 
= [100-62Y"' = 4.2'-^ 
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Time to percolation, p,-p 

Fierure 5: Parameter c of the mass-rank relation function of time. At 

percolation c(pc) ~ 4.2; the Euclidean limit of ^] corresponds to c = 3.25, it is reached 
at pc — P ~ 0.14. 




Mass, m 

Figure 6: Mass distribution of clusters observed on a 2, 000 x 2, 000 lattice at percolation 
p = Pc (dash-dotted line), p = 0.48 (dashed line), and averaged over < p < pc (solid 
line). To smooth out statistical fluctuations we show the cumulative distribution: oc 
J2m'>m''^m- For comparlsou, all curves are normalized to unity at m = 1. 
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Rank, r 



Figure 7: Rank distribution of clusters observed for 2, 000 x 2, 000 lattice at percolation 
p = Pc (dash-dotted line), p = 0.29 (dashed line), and averaged over the percolation cycle 
< p < Pc (solid line). For comparison, all curves are normalized to unity at r = 1. 
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Figure 8: Dynamics of population n,. ■ of a given rank, r = 5, 6, 7 for L = 2, 000. 
Moment of percolation is depicted by a vertical dashed line. 
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0.2 0.4 0.6 0.8 



Density (time), p 

Figure 9: Dynamics of population of a given rank, 1 < r < 11 in semilogarithmic 
scale. Moment of percolation is depicted by a vertical dashed line. (Cf. Fig. 7 from P^). 




Figure 10: Scaling for rank dynamics, a) Ratios Vr^z) j vii^z) do not collapse thus rejecting 
the two-exponent scaling hypothesis; see details in Sect. 14.231 b) Position of percolation 
on the normalized Gaussian qq^z — 0.51); see details in Sect. 14.2.31 
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Figure 11: Three-exponent scaling for rank dynamics. The master Gaussians go{z) for 
different ranks collapse when using the renormalization given by Eqs. ()12|) . (fT^ . 
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Figure 12: Corrections to scaling. The pure exponential rank distribution of Eq. Q 
suggests a horizontal plateau for the normalized function Nr = 10'''" n^, while the observed 
values clearly deviate from the plateau at small and large clusters (panel a). The large 
cluster deviation is due to finite size efffects and is described by an exponential correction 
of Eq. (USD with ^ 0.33 (panel b). 
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Time to percolation, p,-p 



Figure 13: Dynamics of scale corrections. A convexity measure /i(p) is defined by Eq. (pUj) 
and illustrated in the insert. It is positive for convex, and negative for concave rank 
distribution. The downward bend of the rank distribution observed at early stages {fi > 
0) is changed to the upward one (/x < 0) for {p^ — p) < 2 ■ 10~^. See details in Sect. 14.4.21 
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Figure 14: Dynamics of number of clusters rty. (solid line), total mass Mj. (dashed line), 
and average mass (dash-dotted line) for clusters of rank r = 5. 




Figure 15: Fractal structure of clusters, a) Mass-circumference relation for clusters of 
different ranks. The asymptotic power relation with slope 1.0 is gradually develops as 
rank increases, b) Values of fractal dimension D^. (Eq. ((221)) for different ranks. Both 
panels correspond to a 2, 000 x 2, 000 lattice at percolation. 
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Figure 16: Premonitory increase of cluster fractal dimension. The steady-state dynamics 
of fractal dimension (Eq. |221) changes, and starts to increase, as system approaches 
percolation. Similar phenomenon is observed for other ranks. 
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Figure 17: Dynamical constraint for nr{p). Dynamics of 9^ = logln^/n^) and Uq is shown 
in panel a: peaks of two curves coincide. The similar phenomenon is observed for other 
ranks: panel b shows the times of maxima of (circles) and minima of dr-2 (triangles) 
for 3 < r < 9. 
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